function S = consolidate_mmsample(R, thres)
% Consolidate mixture model samples
%
%   S = consolidate_mmsample(R);
%   S = consolidate_mmsample(R, thres);
%
%       summarize the samples of mixture model.
%
%       R is the result of Gibbs sampling (such as those returned by
%       fmm_gibbs, or dpmm_gibbs).
%       
%       S is the consoliated mixture model, which is a struct that 
%       comprises the following fields:
%
%       - 'Id':     the (retained) indices of components (1 x K')
%       - 'Theta':  the matrix of component paremeters (d x K')
%       - 'W':      the (posterior) weights of components (1 x K')
%
%       Here, K' is the number of retained components (which can be
%       less than the total number in the entire system).
%       Theta(:,k) and W(k) corresponds to Ids(k).
%
%       If thres is given, then models whose weights are below thres
%       are discarded.
%

% Create by Dahua Lin, on Nov 16, 2010
%

%% verify input

if ~(isstruct(R) && isfield(R, 'tag') && strcmp(R.tag, 'mm_samples'))
    error('consolidate_mmsample:invalidarg', ...
        'R should be a struct representing mmsample.');
end

if nargin < 2
    thres = 0;
else
    if ~(isfloat(thres) && isscalar(thres) && isreal(thres) && thres >= 0)
        error('consolidate_mmsample:invalidarg', ...
            'thres should be a non-negative real scalar.');
    end
end

%% main

samples = R.samples;

% join

Ids = [samples.Id];
Ws = [samples.W];
Thetas = [samples.Theta];
d = size(Thetas, 1);

% summarize

[Id, gs] = uniqvalues(Ids, 'G');
m = numel(Id);
W = zeros(1, m);
Theta = zeros(d, m);

for k = 1 : m
    gk = gs{k};
    nk = numel(gk);
    W(k) = sum(Ws(gk)) * (1/nk);
    Theta(:,k) = sum(Thetas(:, gk), 2) * (1/nk);
end

% filter (thresholding)

if any(W <= thres)
    si = find(W > thres);
    Id = Id(si);
    W = W(si);
    Theta = Theta(:, si);    
end

S = struct( ...
    'Id', Id, ...
    'Theta', Theta, ...
    'W', W);


